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Phase diagram of QCD at finite temperature and chemical potential from lattice 
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We present the first results for lattice QCD at finite temperature T and chemical potential fi 
with four flavors of Wilson quarks. The calculations are performed using the imaginary chemical 
potential method at K = 0, 0.001, 0.15, 0.165, 0.17 and 0.25, where k is the hopping parameter, 
related to the bare quark mass m and lattice spacing a by « = l/(2ma + 8). Such a method allows 
us to do large scale Monte Carlo simulations at imaginary chemical potential /i = ifii. By analytic 
continuation of the data with /ij < nT/3 to real values of the chemical potential, we expect at each 
k £ [0, Kchiral], a phase transition line on the (fJ,,T) plane, in a region relevant to the search for 
quark gluon plasma in heavy-ion collision experiments. The transition is first order at small or large 
quark mass, and becomes a crossover at intermediate quark mass. 

PACS numbers: 12.38.Gc, ll.10.Wx, 11.15.Ha, 12.38.Mh 



I. INTRODUCTION 

QCD predicts a transition between hadronic matter 
(quark confinement) to quark-gluon plasma (quark de- 
confinement) at sufficient high temperature T and small 
chemical /i. This matter might have existed in early uni- 
verse immediately after the big bang. The main purpose 
of heavy-ion collision experiments at LHC (CERN), SPS, 
and RHIC (BNL) is to recreate such an environment. At 
large \i and lower T, several QCD-inspired models predict 
the existence of a color-superconductivity phase, which 
might be relevant to neutron star or quark star physics. 
Therefore, it is of great significance to study the QCD 
phase structure at larger \i. Because QCD is still strongly 
coupled at criticality, perturbative methods do not apply. 
Lattice gauge theory (LGT) is the most reliable tool for 
investigating the phase transition from first principles. 

In the continuum, the thermodynamics is described by 
the grand partition function 



Z{jjL, T) = Tr e 



-(H-nN q )/T 



(1) 



where H is the Hamiltonian, fi is the chemical potential, 
and N q = J d xtpip is the quark number operator. 

In the Lagrangian formulation of SU(3), LGT at fi- 
nite fi, the effective fermionic action is complex and 
traditional Monte Carlo (MC) techniques with impor- 
tance sampling do not work. The recent years have seen 
enormous efforts H on solving the complex action 
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problem, and some very interesting information^] on 
the phase diagram for QCD with Kogut-Susskind (KS) 
fermions at large T and small fi has been obtained. The 
KS approach to lattice fermions thins the fermionic de- 
grees of freedom in naive fermions, but it does not com- 
pletely solve the species doubling problem. It preserves 
the partial chiral symmetry, but it breaks the flavor sym- 
metry. One staggered flavor corresponds to four flavors 
in reality and the fermionic determinant is replaced by 
the fourth root. Such a replacement is mathematically 
un-provenp^ and it might lead to the locality problem in 
numerical simulations |6|. 

Wilson's approach to lattice fermionsQ avoids the 
species doubling and preserves the flavor symmetry, but 
it explicitly breaks the chiral symmetry, one of the 
most important symmetries of the original theory. Non- 
perturbative fine-tuning of the bare fermion mass has 
to be done, in order to define the chiral limit. There 
have been some MC simulations of LGT with Wilson 
fermions at finite temperature, but no numerical investi- 
gation at finite chemical potential. In Ref. the QCD 
phase structure on the (/i, T) plane was investigated us- 
ing Hamiltonian lattice QCD with Wilson fermions at 
strong coupling. For Nf/N c < 1, with N c = 3 being the 
number of colors and Nf being the number of flavors, a 
tricritical point is found, separating the first order and 
second order chiral phase transitions. 

Neuberger's overlap fermion approach^ solves the 
species doubling problem and preserves the chiral sym- 
metry and flavor symmetry, by introducing exponentially 
decaying non-local terms in the action. However, the 
computational costs[Tol| for dynamical overlap fermions 
are typically two orders of magnitude heavier than for 
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the Wilson or KS formulations. This is certainly beyond 
the current computer capacity. 

In this paper, we perform the first MC simulations of 
Lagrangian LGT with four flavors of dynamical Wilson 
quarks at finite temperature and imaginary chemical po- 
tential an d study the phase structure in the ((a, T, k) 
space, by analytical continuation to the real chemical po- 
tential. 



II. LATTICE FORMULATION 

The basic idea of lattice gauge theory Q, as proposed 
by K. Wilson in 1974, is to replace continuous space- 
time by a discrete grid with lattice spacing a. Gluons 



live on links Uj(x) — e 



/ J 3 dx' Aj (x') 



and quarks live 



on lattice sites. The continuum Yang-Mills action S g 
J d 4 x Tr Fji(x)Fji(x)/2 is replaced by 



(2) 



where (3 = G/g 2 , and U p is the ordered product of 
link variables U around an elementary plaquette. The 
continuum quark action Sf = j d 4 x r ip cont (x)('yjDj + 
m)ip cont (x) is replaced by 



S f = 5^(a;)M B|S V(i/). 



(3) 



For Wilson fermions, the quark field ip on the latt ice is re- 
lated to the continuum one ip cont by ip = ijj cont y / a 3 /(2n) 
with k = 1/ (2ma + 8) the hopping parameter. M is the 
fermionic matrix: 



M^x 7 y — ^x,y 



3=1 



(1 - j 3 )U 3 (x)6 x ^ 



+ (l + -Yj)U]{x-j)6 StU+3 



(4) 



If we consider Nf degenerate flavors of quarks, the par- 
tition function is 



Z = / [dU][dip}[dip}e- s <> 



-s f 



= I [dU] (DetM[U]) Nf e~ s « 



(5) 



In the Lagrangian formulation of LGT, following Ref. 
1 1 1| , the chemical potential is introduced by replacing the 
link variables in the temporal direction in fermion action 
with: 



U 4 (x) ™> e a "U 4 (x), U\{x) -> e- af *Ul(x). 



(6) 



The fermionic action is reduced to the continuum one 
when a — > 0. However, the effective fermionic action 
in the partition function becomes complex, and forbids 



MC simulation with importance samplin g. Several re- 
vised methods, e.g., improved reweighting 12J and imag- 
inary chemical potential [T^ methods, were proposed 
to simulate QCD with KS fermions at finite /i. 

Lattice QCD at imaginary chemical potential does not 
suffer the complex action problem. In this paper, we will 
apply this method to the MC study of the phase structure 
of Wilson fermions. We will measure the expectation of 
the Polyakov loop, chiral condensate and their suscepti- 
bilities. 

The Polyakov loop (Wilson line) is defined as 



P(x) = Tr 



'N t -1 

n 



(7) 



where N t is the number of lattice sites in the temporal 
direction. P(x) is used to scale the interaction strength 
between quarks. In pure gauge theory, non-zero P(x) 
is a signal of quark deconfinement. In practice, we 
measure its expectation value over configurations gener- 
ated in MC simulations with the probability distribution 
(DetM[U]) Nf e- s s/Z. 



The chiral condensate is defined as 



(H>) = 



[dU] [#][#] ipipe 



-S„ — Sr 



1 



ZVN t 



[dU]Tr (M-^U]) (DetM[C/]) JV/ e 



(8) 

where V is the spatial lattice volume. With suitable 
subtraction, it serves as the order parameter of chiral- 
symmetry breaking in the non-perturbatively defined chi- 
ral limit, k = Kchirai where the pion is massless. 

Strictly speaking, if the dynamical quarks play a role, 
P{x) is no longer an order parameter for deconfinement; 
(ipip) is no longer an order parameter for spontaneous 
chiral-symmetry breaking, if k ^ K c hiral- However, when 
the system is at criticality, in particular for the first or- 
der transition, one should observe sharp changes in these 
quantities. Of course, this method is very rough. At the 
transition point, there will be a peak in the susceptibil- 
ity, because the fluctuation of physical quantities is very 
strong. Therefore, the susceptibility will provide more 
useful information about the transition. The susceptibil- 
ity of a quantity 



O = 



VN, 



is defined as 



X = VN t ((0 2 -(0) 2 )). 



(9) 



(10) 



At criticality, the maximum value of x behaves as Xmax oc 
V a , with a the critical exponent. If a — 0, the transition 
is just a crossover; If < a < 1, it is a second order phase 
transition; If a = 1, it is a first order phase transition, 
accompanying the double peak structure in the histogram 
of the quantity O and flip-flops between the two states 
in the MC history. 
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III. PHASE DIAGRAM ON THE (jii,T) PLANE 

Many years ago, Roberge and Weiss (RW)Jj[ made 
the first analytical study of the phase structure of a gauge 
theory with fermions. Replacing [i by i^ij (/ij being a real 
number) and introducing 9 = PjjT, the grand partition 
function has the behavior 

Z(0) = Z(9 + 2n). (11) 

If the gauge group is SU(JV C ), the exact period is 2ir/N c . 
It implies the theory has a Z3 symmetry for N c = 3, 
which is clearly an artifact of imaginary /1 and unphysical. 
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FIG. 1: Schematic phase diagram of QCD on the (/zr,T) 
plane, suggested by Roberge and Weiss. For T > Te, there 
are first order phase transitions at a/ii — 2n(k + 1/2) /N, 
characterized by discontinuity in the Polyakov loop. Here 
N = N c Nt and Nta = 1/T. [There is no phase transition 
across the dashed lines; they just illustrate the location of Te 
or afn = 2%(k + l/2)/N]. 
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first order RW phase transitions between different Z3 sec- 
tors, characterized by the appearance of discontinuity in 
the Polyakov loop. However, later MC simulations with 
KS fermions |0, 0] suggest a phase diagram as shown 
in Fig. |2 i.e., there are additional chiral/deconfincmcnt 
transition lines for Tc < T < Te, with Tc being the 
critical temperature at fi = 0. However, only the critical 
line lies within 9 G [0, n/N c ) and T G [T c , T E ) is relevant 
for the analytical continuation to the real chemical po- 
tential in the physical case. I.e., the imaginary chemical 
potential method works only for fii/T < ir/N c . 

IV. MC SIMULATIONS 

In this section, we will present the first results for QCD 
with four flavors of Wilson fermions at finite T and i/ii. 

The R algorithm^ was used. We modified the MILC 
collaboration's public LGT codej26] to simulate the case 
of imaginary chemical potential. The simulations were 
done at lattice size V x N t — 8 3 x 4 and hopping pa- 
rameter k =0, 0.001, 0.15, 0.165 0.17, and 0.25. At some 
a/i/, j3 and k values, finite size scaling analysis was per- 
formed on different lattices. There are 20 molecular steps 
with size St = 0.02 for each configuration. For each (3 
and afii, we generated at least 20,000 configurations, af- 
ter 4000 warmups. 20 iterations are carried out between 
measurements. Around the area where the thermody- 
namical observables change rapidly, we raised the statis- 
tics at least two times. All simulations were done on 
our PC clusters. The cluster with 20 Pentium III-500 
CPUs[H 13 was built in 2000, and has been upgraded 
to 60 CPUs, with 40 new AMD Opteron-242 CPUs. 
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FIG. 3: Polyakov loop norm as a function of a[ii at k — 0.15 
and two different f3. 



FIG. 2: Schematic phase diagram of QCD on the [fjii,T) 
plane, suggested by lattice MC study. 

Figure^shows the phase diagram of QCD at imagi nary 
chemical potential, suggested by Roberge and Weiss |15|. 
Below some Te, there is no phase transition. Above Te 
and at 6 = 2?r(fc + l/2)/N c , with k = 0, 1,2, there are 



We have more complete data at k = 0.15 than other k 
values. Figure plots the results for the Polyakov loop 
norm (|P(i)|) as a function of a/ij at two different and 
larger values of j3 (corresponding to higher T) . The data 
for (■tpt/j) are shown in Fig. 21 As one sees, these quan- 
tities are approximately periodic with period 7r/6, con- 
firming the RW period for afn to be 27raT/3 = 2n/{3N t ). 
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FIG. 4: Chiral condensate as a function of a^ii at « = 0.15 FIG. 6: MC history of the phase of the Polyakov loop at 
and two different /3. a/ii = 0.26, /3 = 5.25, and k, = 0.15. 



Furthermore, at a/ij = 2(k + l/2)n/(3N t ), k = 0,1,..., 
there is a rapid change in the thermodynamical quan- 
tities, indicating the RW phase transition at higher T 
(larger (3). Figure [3] is a more detailed scan for the phase 
of the Polyakov loop at f3 = 5.25. At a/ij = 0.26, it 
changes rapidly. Figure [S] plots the MC history of the 
phase of the Polyakov loop. There is a clear signal for 
first order phase transition at a/ii = 0.26, confirming the 
existence of RW transition between different Z3 sectors 
at afii = n/(3N t ) = tt/12. 




FIG. 5: Phase of the Polyakov loop as a function of a/i/ at 
P = 5.25 and k = 0.15. 

To locate the chiral/deconfincmcnt phase transition 
line, we made more detailed measurements of (\P(x)\), 
X\p\, (iW, and for a/i/ < w/(3N t ). 

Figures [7| and |H1 show respectively (ifiip) and x$ii> ver ~ 
sus P at a/i/ = 0.05 and 0.20 for k — 0.15; Figures 
and 1101 show the results for (\P\) and X\P\- In the rapid 
changing area, the Polyakov loop and its susceptibility 
behave more singularly than the chiral condensate and 
chiral susceptibility. From the position of the peak in 
the susceptibilities, we determine the transition point. A 
collection of transition points (a/j,i,f3c) is listed in Tab. 
ID In Refs. EEG3, 

it has generally been argued that for 
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FIG. 7: Chiral condensate as a function of f3 at ft : 
two different values of a/ii. 



0.15 and 



small a/ii, the transition line Pc(o-fJ-i) can be expressed 
as a Taylor series with even power of afij . Due to limited 
data, the series is truncated to the quadratic term. We 
use the least squares method to fit the data in Tab. [I]for 
k = 0.15, and obtain an equation for the transition line: 

/3 C = 5.169(9) + 0.954(33) (a/i/) 2 + O (aVf) , (12) 

with error bars coming from the fit. 

Nevertheless, at this k, there is no obvious double peak 
structure in the histograms of thermodynamical observ- 
ables. To determine the nature of the transition, one has 
to do a finite size study. Let us take the transition point 
at a/i] = 0.14 and (3 — 5.187 in Tab. [I]as an example. 
FigurelTTlconrpares x\p\ f° r spatial volumes V = L 3, = 8 3 
and 12 3 around the transition point (with A/3 = 0.002). 
Within error bars, the locations and heights of the peaks 
are consistent. At this f3, we did the longest simulations 
on 8 3 x 4, 10 3 x 4, 12 3 x 4, 14 3 x 4, and 16 3 x 4 lattices, 
with statistics at least four times higher than other non- 
transition points. As shown in Fig. 1121 X\p\ does not 
increase as L. This implies that the transition point at 
a/ij = 0.14 and (3 = 5.187 for k = 0.15 is just a crossover. 
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FIG. 8: Chiral susceptibility as a function of /3 at « = 0.15 
and two different values of a[ii. 




P 

FIG. 9: Polyakov loop norm as a function of /3 at ft = 0.15 
and two different a/xj. 



However, the situations at small or large k are very 
different. Figure IT31 shows the histogram of the Polyakov 
loop norm for (3 = 4.870 and afii = 0.1, at k = 0.165 
which is closer to the chiral limit K c hiral- We observe a 
double peak structure. Figure ITU plots the history of the 
MC simulation at the same parameters, with \P\ mea- 
sured after warmups; We observe the jumps of \P\ from 



a^i Pc 
0.00 5.168(2) 
0.05 5.170(2) 
0.10 5.180(2) 
0.14 5.187(2) 
0.18 5.200(2) 
0.20 5.206(2) 
0.22 5.217(2) 

TABLE I: Collection of transition points for k, — 0.15, deter- 
mined by locating the peak of the susceptibilities, with error 
bars coming from the scan precision. 
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FIG. 10: Susceptibility of the Polyakov loop norm as a func- 
tion of f3 at k = 0.15 and two different afii. 
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FIG. 11: Susceptibility of the Polyakov loop norm as a func- 
tion of P at a/ij = 0.14 and k = 0.15 for different spatial 
volumes V = 8 3 and 12 3 . 



one plateau to another. The results for k = 0, 0.001 and 
0.17 are similar, as shown in Figs. 1151 1161 and 1171 These 
indicate that the phase transitions at k G [0, fti] and 
K G [re 2 , K chirai] are of first order; Here Ki G (0.001, 0.15), 
k 2 G (0.15,0.165) and n chiral G (0.17,0.25). At ft = 0.25, 
which should be the case when k > K c hirah Figs. 1181 and 
EH tell us that there is no phase transition. 



V. PHASE DIAGRAM ON THE (fi,T) PLANE: 
THE PHYSICAL CASE 

Replacing /ij by — iji, we directly continue the transi- 
tion line (|12f> at n = 0.15 from imaginary chemical po- 
tential to real chemical potential: 

P c = 5.169(9) - 0.954(33) (a/i) 2 + O (a 4 fi 4 ) . (13) 

In order to translate the lattice results into the physical 
units, we use the renormalization group relation between 
the lattice spacing a and (3. The two loop perturbative 
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FIG. 12: Susceptibility of the Polyakov loop norm as a func- 
tion of the spacial extent L at a\ii = 0.14, f3 = 5.187 and 
k = 0.15. 




FIG. 13: Histogram of the Polyakov loop norm at and afii 
0.1, = 4.87 and k = 0.165. 



FIG. 14: MC history of the Polyakov loop norm at afii = 0.1, 
(3 = 4.87, and k = 0.165. 
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FIG. 15: MC history of the Polyakov loop norm at afJ,i = 0.1, 
P = 5.691, and k = 0. 



VI. DISCUSSIONS 



expression gives [21 



exp 



4tt 2 



33 - 2N f 



459 - 57iV/ 
(33 - 27V/) 2 



In 



8tt 2 



33 - 2N f 







(14) 
(15) 



where A/, is the lattice QCD scale. In Ref. [22|, the phase 
transition was studied with different kinds of fermion 
actions at fi = and consistent results were obtained. 
Therefore we fix the lattice QCD scale by the critical 
temperature To = 164 MeV at \i = with 4 flavors for 
KS fermions |2l|. The temperature T is related to a and 
N t byT = l/(aN t ). The critical line on the (p,T) plane 
is shown in Fig. I2UI For comparison, the critical line for 
KS fermions|17| is also shown. Within error bars, the 
results are consistent. 



In the preceding sections, we have studied the prop- 
erties of the phase structure of four-flavor QCD on the 
(p, T) plane, using the information obtained from MC 
simulations of LGT with Wilson fermions at imaginary 
chemical potential ifij < iirT/3. 

The advantages of Wilson formulation have been men- 
tioned in the introduction: it is free of species doubling 
and there is one to one correspondence between the fla- 
vors on the lattice and in the continuum. 

Our study suggests that QCD with four flavors of Wil- 
son quarks experiences a first order phase transition from 
the confinement phase to the deconfinement phase at 
small or large quark mass. However, at intermediate 
quark mass, the transition becomes a crossover. The 
properties of the transition are similar to those of KS 
fermions. This region is of interest for the present heavy- 
ion collision experiments. 

Figurel2l"lis the expected phase diagram of lattice QCD 
with Wilson fermions in the (/i, T, k) parameter space. 
There is a surface k = K c hirai where the pion becomes 
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FIG. 16: MC history of the Polyakov loop norm at afj,i = 0.1, FIG. 18: Chiral condensate susceptibility as a function of 
/3 = 5.68, and k = 0.001. at afu = 0.1 and k = 0.25. 





FIG. 17: MC history of the Polyakov loop norm at am =0.1, FIG ' 19: Susceptibility of the Polyakov loop norm as a func- 
P = 4.736, and k = 0.17. tlon of f 3 at a ^ = 0A and K = °- 25 ' 



massless. Above this surface, there is no phase transition, 
as confirmed by our numerical simulations for k = 0.25. 
Interesting physics is below this surface: at each k, one 
should see a phase structure similar to Fig. I^D] Of course, 
the order of transition depends on the value of n. Due to 
heavy computational costs of simulations with dynamical 
fermions, we did not do a comprehensive search for the 
exact location of k±, K2 and K c hiral- We believe that K\ £ 
(0.001,0.15), k 2 6 (0.15,0.165) and n chmd G (0.17,0.25). 

For lower temperature and larger chemical potential, 
all available MC simulation techniques fail. In the Hamil- 
tonian lattice formulation, there has been successful anal- 
ysis of the critical behavior for strong coupling QCD 
lH iU at ((j,, T = 0) and on the (fi,T) planeQ. To 
study the continuum physics, new methods have to be 
developed. We hope to study these issues, as well as the 
dependence of the phase structure on the quark flavors, 
in the near future. 
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FIG. 20: Phase diagram on the {n,T) plane for k = 0.15. 
The area between the two solid lines is our result for Wilson 
quarks with error band, derived from Eqs. I|13|l and 114|l . The 
area between the two dotted lines is the result for KS fermions 
by D'Elia and Lombardo. 
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FIG. 21: Expected phase diagram of lattice QCD with four 
flavors of Wilson quarks in the (fj,, T, k) parameter space. For 
k £ [0, Ki] and k £ [«2, Kchiral], the phase transition is of first 
order, and while for k € K2), the transition is a crossover. 
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